function [fig_handle, h] = gmmplot(neuralGasMap, mmu, ssigma, plotDist, fig_handle);

% GMMPLOT Plotting of a Gaussian Mixture contour
%
%    [h] = gmmplot(X, mmu, ssigma, fig_handle);
%    GMMPLOT (MU,SIGMA) plots the mean and standard
%    deviation ellipsis of the Gaussian processes with
%    means mmu and covariances ssigma
%    
%    mmu  (1 x k) cell array of (1 x dim) mean vectors
%         OR
%         (NUM_KERNELS x dim) array of means 
%    ssigma (1 x k) cell array of (dim x dim) cov matrices
%           OR
%           (NUM_KERNELS * dim^2) array of cholesky factorizations
%
%    further assumes:
%    global COLORS   (NUM_KERNELS+1 * 3)   global color definition
  
  global COLORS

  X = neuralGasMap.D;
  
  if nargin < 5
    fig_handle = figure;
    %plot(X(:,1),X(:,2),'+','color',[.5 .5 .5]), hold on;
  else
    % assume hold on and data already plotted
    fig_handle = figure(fig_handle);
  end

  npts = 100;
  ncircles = 7;
  
  t = linspace(-pi, pi, npts);
  t=t(:);
  
  if iscell(mmu)
    mmu_c = mmu;
    clear mmu;
  
    mmu = cat(1,mmu_c{:});
  end
  
  [NUM_KERNELS dim] = size(mmu);

  for k = 1:NUM_KERNELS
    if iscell(ssigma)
      sigma = ssigma{k};
    else
      R = reshape(ssigma(k,:),dim,dim);
      sigma = R'*R;
    end
    
    stdev(:,:,k) = sqrtm(sigma);
  end
  

  switch dim
   case 2
    plot(X(:,1),X(:,2),'x','color',[1 0 0]), hold on;
    if ~isempty(neuralGasMap) 
        M=neuralGasMap.M;
        linkMatrix=neuralGasMap.linkMatrix;
        plot(M(:,1),M(:,2),'.','color',[0 0 0]), hold on;
        
        %linkmatrix
        bmus=1:size(M,1);
        for i=1:size(linkMatrix,1)
            links = linkMatrix(i,:);

            linkIndex = links >= 0;
            linkedBmus=bmus(linkIndex);
            for j=1:length(linkedBmus)
                 bmu1=M(i,:);
                 bmu2=M(linkedBmus(j),:);
                 x=[bmu1(1) bmu2(1)];
                 y=[bmu1(2) bmu2(2)];
                 line(x,y);
            end;
        end;
    end;
    %pause;   
    hold off
    
    if (plotDist == 1)
        for k = 1:NUM_KERNELS
          SIG = [cos(t) sin(t)] * stdev(:,:,k) + repmat(mmu(k,:),npts,1);

          h(k,1) = line(SIG(:,1),SIG(:,2),'color',COLORS(k,:),'linew',2);
          h(k,2) = line(mmu(k,1),mmu(k,2),'marker','+','markersize',10,'color',COLORS(k,:),'linew',2); 
        end
    end
    % stay with hold on if figure was provided
    %if nargin < 4, hold off; end
     
    hold off, drawnow;

   case 3
    plot3(X(:,1),X(:,2),X(:,3),'+','color',[.7 .7 .7]), hold on;
    
    % see plotgaus3 for details
    coeff = linspace(0.9, -0.9, ncircles)';
    coeff(:,2) = sqrt(1-abs(coeff(:,1)).^2);
    X = [];

    for j = 1:ncircles
      X = [X; [coeff(j,2)*cos(t) coeff(j,2)*sin(t) coeff(j,1)*ones(npts,1)]];
    end

    for k = 1:NUM_KERNELS
      X_ = X*stdev(:,:,k);
      X_ = X_ + repmat(mmu(k,:), npts*ncircles, 1);

      h(k,1) = plot3(X_(:,1), X_(:,2), X_(:,3), '.','color',COLORS(k,:),'linew', 2);
      h(k,2) = line(mmu(k,1),mmu(k,2),mmu(k,3),'marker','+','markersize',20,'color', [0 0 0],'linew',1);
    end
    
    hold off, drawnow;
    %fprintf('\n press any key to continue\n'); pause;
    
   otherwise
    fprintf('\n Only 2 or 3 dimensional data supported \n');
  end
  